rm(list=ls())
library(metafor)
library(meta)
library(dmetar)
library(puniform)
library(ivreg)
library(ivpack)
###################
#### Loading Data and AK Estimators
###################
PATH="C:/Users/hong_/Desktop/World Development/FinalCodes_Data/"
DT<-as.data.frame(read.csv(paste(PATH, "data_sorted.csv", sep='')))


### P-uniform Test (publication bias exists i)
# L.pb = test statistic of p-uniforms publication bias test
# pval.pb = one-tailed p-value of p-uniforms publication bias test
puniform(ni = DT$df, tobs = DT$pcc/DT$sepcc, method="LNP",
         side = "right", plot = FALSE)


### Instrument regressions
pubbias_iv = ivreg(pcc ~ sepcc | df, weights=(1/sepcc^2), data=DT)
summary(pubbias_iv)
cluster.robust.se(pubbias_iv, DT$id)


### P-Curve (requires R version higher than 3.5.2)
#install.packages("devtools")
#devtools::install_github("MathiasHarrer/dmetar")
library(dmetar)
meta_obj <- metagen(DT$pcc,DT$sepcc)
pcurve(meta_obj)



### Caliper Test
## Comparing frequencies of p-values in an interval just below .05 
## to the frequency in the adjacent lower interval; 
## a higher frequency in the interval closest to .05 
## is evidence for QRPs that seek to obtain just significant results. 
DT$abtstat <- abs(DT$pcc/DT$sepcc)
range <- c(sort(c(0,seq(1.96,0, -0.20))),sort(seq(1.96+0.2,90, 0.20)))
col <- rep('white', length(range))
col[c(10,11)] <- 'orange'
x11()
hist(DT$abtstat,
    breaks=c(sort(c(0,seq(1.96,0, -0.20))),sort(seq(1.96+0.2,90, 0.20))),
    xlim=c(0.16,13.16),
    freq=TRUE,
    main='Histogram of t-statistics',
    xlab='abs(t-statistics)',
    col=col
)
abline(v=1.96, lwd=2, col='red')
box()

CaliperTable <- as.data.frame(matrix(NA, nrow=3,ncol=4))
colnames(CaliperTable) <- c('Test','OverCaliper','UnderCaliper','p-value')
CaliperTable[1, ] <- c('10% Caliper',sum(1.96<DT$abtstat & DT$abtstat<2.16), sum(1.76<DT$abtstat & DT$abtstat<1.96), NA)
CaliperTable[2, ] <- c('15% Caliper',sum(1.96<DT$abtstat & DT$abtstat<2.26), sum(1.66<DT$abtstat & DT$abtstat<1.96), NA)
CaliperTable[3, ] <- c('20% Caliper',sum(1.96<DT$abtstat & DT$abtstat<2.36), sum(1.56<DT$abtstat & DT$abtstat<1.96), NA)
CaliperTable[, 4] <- c(binom.test(as.numeric(c(CaliperTable$OverCaliper[1],CaliperTable$UnderCaliper[1])), alternative = c("greater"))$p.value,
                        binom.test(as.numeric(c(CaliperTable$OverCaliper[2],CaliperTable$UnderCaliper[2])), alternative = c("greater"))$p.value,
                        binom.test(as.numeric(c(CaliperTable$OverCaliper[3],CaliperTable$UnderCaliper[3])), alternative = c("greater"))$p.value)
print(CaliperTable)







